#ifndef PREPROCESS_CL
#define PREPROCESS_CL
#include <clUtil.cl>
#include <tables.cl>
#include <constants.cl>

__kernel void preprocessSpectra(__global unsigned int* intensity,
                                __global unsigned int* mz,
                                __global unsigned int* intensity2,
                                __global unsigned int* mz2,
                                unsigned int stride,
                                __global MSScanInfo* scanInfo)
{
  __local unsigned int scratch1;
  __local unsigned int scratch2;
  __local unsigned int scanZeros[64];
  __local unsigned int scanOnes[64];
  size_t mySpectra = get_global_id(1);
  size_t threadID = get_local_id(0);
  unsigned int arrayLen;
  bool isGoodSpectra;

  //Don't bother with full scans
  if(scanInfo[mySpectra].mMSLevel != 1)
  {

    //Update pointers to work on our work item's data
    intensity2 = &intensity2[mySpectra * stride];
    intensity = &intensity[mySpectra * stride];
    mz = &mz[mySpectra * stride];
    mz2 = &mz2[mySpectra * stride];

    arrayLen = scanInfo[mySpectra].mNumPeaks;

    if(arrayLen == 0)
    {
      scanInfo[mySpectra].mIsGoodSpectra = false;
      return;
    }

    arrayMax(intensity, &scratch1, arrayLen);

    //Compute total ion current
    sum(intensity, &scratch1, arrayLen);
    barrier(CLK_LOCAL_MEM_FENCE);
    scanInfo[mySpectra].mTotalCurrent = scratch1;
      
    if(threadID == 0)
    {
      unsigned int maxMZ = mz[arrayLen - 1];
      unsigned int minMZ = mz[0];

      scanInfo[mySpectra].mMinMZ = minMZ;
      scanInfo[mySpectra].mMaxMZ = maxMZ;
      scanInfo[mySpectra].mNumHoles = (maxMZ - minMZ) / 
        (2000.0f * kFragmentMZTolerance);
    }

    //Sort greatest to least by intensity
    radixSortGL(intensity, 
                mz, 
                intensity2,
                mz2,
                arrayLen, 
                &scratch1, 
                &scratch2,
                scanZeros,
                scanOnes);

    if(threadID == 0)
    {
      float cutoff = 
        (float)kTICCutoffFactor * (float)scanInfo[mySpectra].mTotalCurrent;
      float currentSum = 0.0f;

      //Find cutoff index
      for(unsigned int curPeak = 0; curPeak < arrayLen; curPeak++)
      {
        currentSum += intensity2[curPeak];

        if(currentSum > cutoff)
        {
          scratch1 = curPeak;
          scanInfo[mySpectra].mNumPeaks = curPeak;
          break;
        }
      }

      //Make sure we don't exceed the max peak count parameter
      scratch1 = min(scratch1, kMaxPeakCount);

      //Round down the number of peaks to fit in class model
      scratch1 = (scratch1 / ((1 << kNumClasses) - 1)) 
        * ((1 << kNumClasses) - 1);
    }


    //Everybody wait on the guy computing the cutoff...
    barrier(CLK_LOCAL_MEM_FENCE);

    //Everybody fetch the new truncated array size
    arrayLen = scratch1;

    //If this spectra is shit, mark it as such
    isGoodSpectra =  arrayLen < (1 << kNumClasses) - 1 ? false : true;

    if(threadID == 0)
    {
      scanInfo[mySpectra].mIsGoodSpectra = isGoodSpectra;
    }

    if(isGoodSpectra == true)
    {
      //Compute class threshholds where each class is 2x the previous
      if(threadID == 0)
      {
        //2^numClasses - 1
        unsigned int sumBinWeights = (1 << kNumClasses) - 1;
        unsigned int binSize = arrayLen / sumBinWeights;
        unsigned int curIdx = 0;
        unsigned int curPeak = 0;

        for(unsigned int curClass = 0; curClass < kNumClasses; curClass++)
        {
          curIdx += (1<<curClass) * binSize;

          //Overwrite the intensity field with its class
          for(; curPeak < curIdx; curPeak++)
          {
            intensity2[curPeak] = curClass;
          }
        }
      }

      //Sort least to greatest by m/z
      radixSortLG(mz2,
                  intensity2,
                  mz,
                  intensity,
                  arrayLen,
                  &scratch1,
                  &scratch2,
                  scanZeros,
                  scanOnes);

      //Fill in minMZ, maxMZ, and num holes entries
      if(threadID == 0)
      {
        scanInfo[mySpectra].mNumPeaks = arrayLen;
      }
    }
  }
}

#endif

